#assembling_binning_features
system("wget ftp://genome-ftp.stanford.edu/pub/yeast/protein_info/protein_properties.tab")
proteinProperties <-
    read.table("protein_properties.tab",
               header = FALSE, sep = "\t")
head(proteinProperties)
nms <-
    readLines(textConnection("FEATURE (ORF name)
SGDID
MOLECULAR WEIGHT (in Daltons)
PI
CAI (Codon Adaptation Index)
PROTEIN LENGTH
N TERM SEQ
C TERM SEQ
CODON BIAS
ALA
ARG
ASN
ASP
CYS
GLN
GLU
GLY
HIS
ILE
LEU
LYS
MET
PHE
PRO
SER
THR
TRP
TYR
VAL
FOP SCORE (Frequency of Optimal Codons)
GRAVY SCORE (Hydropathicity of Protein)
AROMATICITY SCORE (Frequency of aromatic amino acids: Phe, Tyr, Trp)
Feature type (ORF classification: Verified, Uncharacterized, Dubious)"))
nms
longnames <- nms
nms[c(1:9,30:33)] <-
    c("yORF",
      "SGDID",
      "molwt",
      "pi",
      "cai",
      "length",
      "nterm",
      "cterm",
      "codonBias",
      "fop",
      "gravy",
      "aromaticity",
      "type")
names(proteinProperties) <- nms
head(proteinProperties)
rownames(proteinProperties) <- proteinProperties$yORF
save(proteinProperties, file = "proteinProperties.rda", compress = TRUE)
promptData(proteinProperties)
head(proteinProperties)


#save(proteinProperties, file = "proteinProperties.rda", compress = TRUE)

features <- subset(proteinProperties,
  select= c("molwt", "pi", "cai", "length",
    "codonBias", "fop", "gravy", "aromaticity"))


# explorative analysis, looking for suitable variable transformation

print(summary(features))

# keeping the complete cases
nomiss <- complete.cases(features)
print(table(nomiss))
i <- which(nomiss == TRUE)
features <- features[i, ]
cat(str(features))

attach(features)
transformed <- data.frame(log(molwt), log(length),
  sqrt(aromaticity), log(cai))
summary(transformed)

features <- data.frame(features, transformed)

# exploratory data analysis pre-binning
do.hist <- function() {
    par(mfrow=c(2,2))
    for (i in 1:ncol(features)) {
      hist(features[,i], breaks="scott", prob=TRUE,
        main="", xlab=names(features)[i])
      m = mean(features[,i])
      s = sd(features[,i])
      curve(dnorm(x, m, s), col=2, add=TRUE)
      qqnorm(features[,i], main=names(features)[i], cex=.5)
      qqline(features[,i], col=2)
    }
    par(mfrow=c(1,1))
}

#do.hist()  #comment if not needed

# computes the breaks using Scott's rule
# bins the data using cut
#


calc.breaks <- function(x, eps=1e-7) {
  #comment print statements if not needed
  print(match.call())
  if (!is.vector(x)) return (NA)
  n <- nclass.scott(x)
  m <- min(x) - eps
  M <- max(x) + eps
  b <- seq(m, M, length=n+1)
  print(b)
  b
}


# binning steps - calculate breaks using Scott's rule
all.breaks <- list(
    b.molwt = calc.breaks(log(molwt)),
    b.length = calc.breaks(log(length)),
    b.cai = calc.breaks(log(cai)),
    b.aro = calc.breaks(sqrt(aromaticity)),
    b.fop = calc.breaks(fop),
    b.cB = calc.breaks(codonBias),
    b.pi = calc.breaks(pi),
    b.gravy = calc.breaks(gravy))


MRVersion.bin.counts <- function(x, breaks, labels=NULL, ...) {
  #comment print statements if not needed
  print(match.call())
  nclass <- length(breaks)-1
  ifelse (is.null(labels),
    labs <- paste("bin", 1:nclass, sep=""), labs <- labels)
  tab <- table(cut(x, breaks=breaks, labels=labs, ...))
  print(tab)
  return(tab)
}


bin.counts <- function(x, breaks, labels=NULL, ...) {
  #comment print statements if not needed
  print(match.call())
  nclass <- length(breaks)-1
  ifelse (is.null(labels),
    labs <- paste("bin", 1:nclass, sep=""), labs <- labels)
  tab <- table(cut(x, breaks=breaks, labels=labs, ...))
  intervals<-as.vector(unlist(labels(table(cut(x, breaks=breaks, ...)))))
  #as.vector(unlist(labels(table(cut(log(molwt), breaks=all.breaks$b.molwt)))))
  print(tab)
  print(intervals)
  return(tab)
  return(intervals)
}


# binning step - bin data
all.counts <- list(
    c.molwt = bin.counts(log(molwt), breaks=all.breaks$b.molwt),
    c.length = bin.counts(log(length), breaks=all.breaks$b.length),
    c.cai = bin.counts(cai, breaks=exp(all.breaks$b.cai)),
    c.aro = bin.counts(aromaticity, breaks=(all.breaks$b.aro)^2),
    c.fop = bin.counts(fop, breaks=all.breaks$b.fop),
    c.cB = bin.counts(codonBias, breaks=all.breaks$b.cB),
    c.pi = bin.counts(pi, breaks=all.breaks$b.pi),
    c.gravy = bin.counts(gravy, breaks=all.breaks$b.gravy))


bin.generate <- function(x, breaks, ...) {
  #comment print statements if not needed
  print(match.call())
  nclass <- length(breaks)-1

  #tab <- table(cut(x, breaks=breaks, ...))
  intervals<-as.vector(unlist(labels(table(cut(x, breaks=breaks, ...)))))

  return(intervals)
}


all.Intervals <- list(
    c.molwt = bin.generate(log(molwt), breaks=all.breaks$b.molwt),
    c.length = bin.generate(log(length), breaks=all.breaks$b.length),
    c.cai = bin.generate(cai, breaks=exp(all.breaks$b.cai)),
    c.aro = bin.generate(aromaticity, breaks=(all.breaks$b.aro)^2),
    c.fop = bin.generate(fop, breaks=all.breaks$b.fop),
    c.cB = bin.generate(codonBias, breaks=all.breaks$b.cB),
    c.pi = bin.generate(pi, breaks=all.breaks$b.pi),
    c.gravy = bin.generate(gravy, breaks=all.breaks$b.gravy))

#createLUT -preparation 
Features<-features[,c(2,3,5,6,7,8,9,10)]
names(Features)[7]<-"molwt"  # consider to do something like if(grep('log.',nam)){nam<-sub('log.', '', nam, perl = TRUE)}
names(Features)[8]<-"length"





#createLUT
#questa funzione va generalizzata per permettere altri names per le features.....non solo pi, gravy codonBias etc....
# df is Features, variants is all.breaks
createLookUpTableForFeatures<- function(featuresdf,variants){

if (!is.list(variants)) return (NA)
if (!is.data.frame(df)) return (NA)

df.cuts <- list(
    c.molwt = cut(df[,"molwt"], breaks=variants$b.molwt),
    c.length = cut(df[,"length"], breaks=variants$b.length),
    c.cai = cut(df[,"cai"], breaks=exp(variants$b.cai)),
    c.aromaticity = cut(df[,"aromaticity"], breaks=(variants$b.aro)^2),
    c.fop = cut(df[,"fop"], breaks=variants$b.fop),
    c.codonBias = cut(df[,"codonBias"], breaks=variants$b.cB),
    c.pi = cut(df[,"pi"], breaks=variants$b.pi),
    c.gravy = cut(df[,"gravy"], breaks=variants$b.gravy))



mydf = list()
envFeat= list()
myLookupF= list()
for(i in 1:length(names(df)))
{
nam<-names(df)[i]
Var2<-paste("c",nam,sep=".")
cat("questa e la feat: ","\n",nam,"\n","questi sono i breaks: ","\n",Var2,"\n")
intervalsAskeys = levels(df.cuts[[Var2]])
moreLetters<-paste(LETTERS, 1:length(LETTERS), sep=".")
values = c(letters,LETTERS,moreLetters)
values<-values[1:length(intervalsAskeys)]
cat("questi sono i values: ","\n",values,"\n")
mydf[[nam]] = data.frame (intervalsAskeys, values, stringsAsFactors=FALSE)
head(mydf[[nam]])
myLookupF[[nam]] = function (key) d [match (key, d [,1]), 2]
envFeat[[nam]] = new.env()
environment (myLookupF[[nam]]) = envFeat[[nam]]
envFeat[[nam]]$d = mydf[[nam]]

}

}

createLookUpTableForFeatures(Features,all.breaks)



#assignAlleles - preparation
#Individuals, Loci,

#assignAlleles
#this function takes as input 2 data frames and one list, the first data frame representing objects #the second representing features, associated to these objects, the list representing the results of #binning the features values; using a methaphor from genetics, objects can be seen as diploid #individuals, features can be seen as the loci of the individual diploid genomes; individuals have a #particular set of allelic variants in the loci of their genomes, constituting their genetic background;
#going more into details of this particular setting for protein interaction prediction, object or diploid #individual here is a protein pair and consequentely object or diploid individual is a putative relation #between two proteins.

#before I need to define this function to get the pair of bins associate to each protein pairs
merge.2columnsIn1 <- function(x,y,z) {

new.column<-numeric(dim(x)[1])

for(i in 1:dim(x)[1]){
new.column[i]<-paste(sort(c(as.vector(x[i,y]),as.vector(x[i,z]))),collapse="")
cat(new.column[i],"\n")
}

return(new.column)
}
#

assignAlleles <- function(individuals, loci, variants, LUT) {
  #comment print statements if not needed
  print(match.call())
  if (!is.list(variants)) return (NA)
  if (!is.data.frame(individuals)) return (NA)
  if (!is.data.frame(loci)) return (NA)
  if (!is.environment(LUT)) return (NA)

colnames(individuals)[1]<-yORF1
colnames(individuals)[2]<-yORF2

Ids<-1:length(loci$yORF)
names(Ids)<-loci$yORF

from<-individuals$yORF1
to<-individuals$yORF2

length(from)
length(to)
length(Ids)
l.from<-length(from)
l.to<-length(to)

vec.from<-Ids[as.character(from)] #Ids has been defined in the pre-binning.R script
vec.to<-Ids[as.character(to)]
from.features<-Features[vec.from,]# nb Features not features
to.features<-Features[vec.to,]
from.features.cuts <- list(
    c.molwt = cut(from.features[,"molwt"], breaks=variants$b.molwt),
    c.length = cut(from.features[,"length"], breaks=variants$b.length),
    c.cai = cut(from.features[,"cai"], breaks=exp(variants$b.cai)),
    c.aromaticity = cut(from.features[,"aromaticity"], breaks=(variants$b.aromaticity)^2),
    c.fop = cut(from.features[,"fop"], breaks=variants$b.fop),
    c.codonBias = cut(from.features[,"codonBias"], breaks=variants$b.codonBias),
    c.pi = cut(from.features[,"pi"], breaks=variants$b.pi),
    c.gravy = cut(from.features[,"gravy"], breaks=variants$b.gravy))
   
#from.features.cuts[["c.cai"]][1:15]
to.features.cuts <- list(
    c.molwt = cut(to.features[,"molwt"], breaks=variants$b.molwt),
    c.length = cut(to.features[,"length"], breaks=variants$b.length),
    c.cai = cut(to.features[,"cai"], breaks=exp(variants$b.cai)),
    c.aromaticity = cut(to.features[,"aromaticity"], breaks=(variants$b.aromaticity)^2),
    c.fop = cut(to.features[,"fop"], breaks=variants$b.fop),
    c.codonBias = cut(to.features[,"codonBias"], breaks=variants$b.codonBias),
    c.pi = cut(to.features[,"pi"], breaks=variants$b.pi),
    c.gravy = cut(to.features[,"gravy"], breaks=variants$b.gravy))

vec.labels.feat=list()
vec2.labels.feat=list()
temp.feat.df=list()
for(i in 1:length(from.names))
{

nam<-from.names[i]
#var<-paste("b",nam,sep=".")
var2<-paste("c",nam,sep=".")

#
vec.labels.feat[[i]]=numeric(l.from)#46 (5000)#1816)#numeric(1135)
for(j in 1:length(from.features.cuts[[var2]])){
cat("questa e int: ","\n",as.character(from.features.cuts[[var2]][j]),"\n")
vec.labels.feat[[i]][j]<-myLookupF[[i]](as.character(from.features.cuts[[var2]][j]))
}
#
vec2.labels.feat[[i]]=numeric(l.to)#46 (5000)#1816)#numeric(1135)
for(j in 1:length(to.features.cuts[[var2]])){
cat("questa e int: ","\n",as.character(to.features.cuts[[var2]][j]),"\n")
vec2.labels.feat[[i]][j]<-myLookupF[[i]](as.character(to.features.cuts[[var2]][j]))
}
temp.feat.df[[i]]<-data.frame(row.names(from.features),from.features[,nam],unlist(vec.labels.feat[[i]]),row.names(to.features),to.features[,nam],unlist(vec2.labels.feat[[i]])) 
}
col.merged=list()
events.df=list()
for(i in 1:length(temp.feat.df)){
 col.merged[[i]] = merge.2columnsIn1(temp.feat.df[[i]],3,6)
 events.df[[i]]<-data.frame(temp.feat.df[[i]],unlist(col.merged[[i]]))
}
for(i in 1:length(events.df))
{
cat("questa e il ciclo: ","\n",i,"\n")
nam<-names(from.features)[i]
cat("this is the feature: ","\n",nam,"\n")
names(events.df[[i]])[7]<-as.character(nam)
}

all.df<-data.frame(events.df[[1]],events.df[[2]],events.df[[3]],events.df[[4]],events.df[[5]],events.df[[6]],events.df[[7]],events.df[[8]])
dim(all.df)

}
#loci <- proteinProperties
dream.df<-read.table("datasets/GoldStandard_ProtProtSubnet.txt")
assignAlleles(dream.df, Features, all.breaks, myLookupF)
